/*
 * Copyright (c) 2007 Thomas Fors.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

// =======================================================
// Variable-size N-dimensional simplex optimization.
// =======================================================



// -------------------------------------------------------
// Simplex()
// Constructor
// -------------------------------------------------------
function Simplex(nDim, name) {
	var i;
	
	this.name = name;
	this.nDim = nDim;					// number of dimensions
	
	this.uv = new Array();				// array of vertices in user coordinates
	this.v = new Array();				// array of vertices in normalized coordinates
	this.r = new Array();				// array of results
	
	this.rCache = new Array();			// result cache
	this.rCacheSize = 0;				// cache size
	this.rCacheHits = 0;				// number of cache hits
	
	this.bestSettings = new Array();	// best settings found
	
	// Construct vertex arrays
	for ( i=0; i<this.nDim+1+4; i++ ) {
		this.uv[i] = new Array(this.nDim);
		this.v[i] = new Array(this.nDim);
		this.r[i] = new Array(this.nDim);
	}
	
	i = this.nDim+1;
	this.R = i++;				// index of Reflected vertex
	this.E = i++;				// index of Expanded vertex
	this.Cr = i++;				// index of Contracted vertex on Wastebasket side
	this.Cw = i++;				// index of Contracted vertex on Reflected side
	
	this.W = -1;				// index of Wastebasket vertex
	this.N = -1;				// index of Next-worst vertex
	this.B = -1;				// index of Best vertex
	
	this.uLo = new Array();		// low boundary in user-space for each dimension
	this.uHi = new Array();		// high boundary in user-space for each dimension
	this.uRes = new Array();	// resolution required in user-space for each dimension
	this.fnSet = new Array();	// array of functions to set each parameter (dimension)
	
	// Load cache
	if ( display.files["cache"] != undefined ) {
		var cacheIndex, cacheValue;
		display.files["cache"].open("e");
		while ( !display.files["cache"].eof ) {
			cacheIndex = display.files["cache"].readln();
			cacheValue = display.files["cache"].readln();
			if ( cacheIndex.length > 0 ) {
				this.rCache[cacheIndex] = cacheValue;
				this.rCacheSize++;
			}
		}
		display.files["cache"].close();
	}
}


// -------------------------------------------------------
// SetDimension()
// Sets user-space boundary and resolution requirements
// -------------------------------------------------------
Simplex.prototype['SetDimension'] = function (d, uLo, uHi, resolution, fnSet) {
	this.uLo[d] = uLo;
	this.uHi[d] = uHi;
	this.uRes[d] = resolution;
	this.fnSet[d] = fnSet;
}

// -------------------------------------------------------
// NormalizeVertex()
// Convert vertex from user space to normalized space
// -------------------------------------------------------
Simplex.prototype['NormalizeVertex'] = function (i) {
	var d, m, b;
	for ( d=0; d<this.nDim; d++ ) {
		m = 1/(this.uHi[d]-this.uLo[d]);
		b = m * (0 - this.uLo[d]);
		this.v[i][d] = this.uv[i][d] * m + b;
	}
}

// -------------------------------------------------------
// DenormalizeVertex()
// Convert vertex from normalized space to user space
// -------------------------------------------------------
Simplex.prototype['DenormalizeVertex'] = function (i) {
	var d, m, b;
	for ( d=0; d<this.nDim; d++ ) {
		m = this.uHi[d]-this.uLo[d];
		b = this.uLo[d];
		this.uv[i][d] = Math.round( (this.v[i][d] * m + b) / this.uRes[d] ) * this.uRes[d];
	}
}

// -------------------------------------------------------
// ProjectVertex()
// Project vertex (i) of simplex through average of remaining
// vertices by amount (factor).
// -------------------------------------------------------
Simplex.prototype['ProjectVertex'] = function (i, factor) {
	var d, j, center;
	var newVertex = new Array();
	
	// for each dimension...
	for ( d=0; d<this.nDim; d++ ) {
		// compute center of remaining vertices
		center = -this.v[i][d];
		for ( j=0; j<this.nDim+1; j++ ) {
			center += this.v[j][d];
		}
		center /= this.nDim;	// num of remaining vertices
		
		// compute new vertex
		newVertex[d] = (1-factor)*center + (factor)*this.v[i][d];
	}
	
	return newVertex;
}

// -------------------------------------------------------
// ComputeReflections()
// Computes R, E, Cr, and Cw vertices from W
// -------------------------------------------------------
Simplex.prototype['ComputeReflections'] = function () {
	this.v[this.R]  = this.ProjectVertex(this.W, -1);
	this.v[this.E]  = this.ProjectVertex(this.W, -2);
	this.v[this.Cr] = this.ProjectVertex(this.W, -0.5);
	this.v[this.Cw] = this.ProjectVertex(this.W, 0.5);
}

// -------------------------------------------------------
// RankVertices()
// Ranks the vertices assigning W, N, and B
// Larger numbers (more positive) are better
// -------------------------------------------------------
Simplex.prototype['RankVertices'] = function () {
	var i;
	var worst;
	var nextWorst;
	var best;
	
	// Initialize W, N, and B before searching
	best = 0;
	if ( this.r[0] > this.r[1] ) {
		worst = 0;
		nextWorst = 1;
	} else {
		worst = 1;
		nextWorst = 0;
	}
	
	// Search...
	for ( i=0; i<this.nDim+1; i++ ) {
		if ( this.r[i] < this.r[worst] ) {
			nextWorst = worst;		// new Next-worst
			worst = i;				// new Worst
		} else if ( this.r[i] < this.r[nextWorst] && i != worst ) {
			nextWorst = i;			// new Worst-next
		}
		if ( this.r[i] >= this.r[best] ) {
			best = i;				// new Best
		}
	}
	
	// Make the assignments
	if ( this.W < 0 ) {
		// we're ranking vertices of the first simplex, so we assign W
		// to the absolute worst point
		this.W = worst;
		this.N = nextWorst;
	} else {
		// otherwise W will have been assigned to the previous N in
		// Step() so we honor that here and just assign N (and B).
		if ( this.W != worst ) {
			this.N = worst;
		}	else {
			this.N = nextWorst;
		}
	}
	// B is always best
	this.B = best;
	
	// Check for new best result and save
	if ( this.r[this.B] > this.bestResult ) {
		this.bestResult = this.r[this.B];
		this.DenormalizeVertex(this.B);
		for ( d=0; d<this.nDim; d++ ) {
			this.bestSettings[d] = this.uv[this.B][d];
		}
	}
}

// -------------------------------------------------------
// DefineOptimization()
// type = -1 for minmization,
//   0 for zero search (untested),
//  +1 for maximization (untested).
// fnEval must accept an array of nDim control values and return a result
// Optimization will return when best-case result beats goal
//   or timeLimit (in minutes) has been exceeded
//   or stepLimit has been exceeded
//   or restartLimit has been exceeded
// if any of these limits are set to undefined then they are disabled
// -------------------------------------------------------
Simplex.prototype['DefineOptimization'] = function (type, fnEval, goal, timeLimit, stepLimit, restartLimit) {
	this.optType = type;
	this.fnEval = fnEval;
	this.goal = goal;
	this.timeLimit = ( timeLimit == undefined ) ? (Number.POSITIVE_INFINITY) : (timeLimit);
	this.stepLimit = ( stepLimit == undefined ) ? (Number.POSITIVE_INFINITY) : (stepLimit);
	this.restartLimit = ( restartLimit == undefined ) ? (Number.POSITIVE_INFINITY) : (restartLimit);
}

// -------------------------------------------------------
// EvaluateVertex()
// Evaluate vertex and record result.  If vertex is out
// of bounds, return highly undesirable (Number.NEGATIVE_INFINITY)
// score.
// -------------------------------------------------------
Simplex.prototype['EvaluateVertex'] = function (i) {
	var d, result;
	var cacheIndex;
	
	//  Test for out of bounds condition
	var bOutOfBounds = 0;
	for ( d=0; d<this.nDim; d++ ) {
		if ( this.v[i][d] < 0 || this.v[i][d] > 1 ) {
			bOutOfBounds = 1;
		}
	}
	
	this.DenormalizeVertex(i);
	
	var bShowStatus = true;
	if ( bOutOfBounds ) {
		this.r[i] = Number.NEGATIVE_INFINITY;	// return extremely bad result to force a contraction
		bShowStatus = false;
	} else {
		// check cache first.  if not present, do the more time-consuming evaluation
		cacheIndex = "";
		for ( d=0; d<this.nDim; d++ ) {
			cacheIndex += this.uv[i][d] + " ";
		}
		result = this.rCache[cacheIndex];
		if ( result == undefined ) {
			for ( d=0; d<this.nDim; d++ ) {
				acr[this.fnSet[d]](this.uv[i][d]);	// make settings
			}
			result = cc[this.fnEval]();				// evaluate
			this.rCache[cacheIndex] = result;
			this.rCacheSize++;
		} else {
			this.rCacheHits++;
			bShowStatus = false;
		}
		
		// count vertex evaluations for degenerate detection
		try {
			(this.evalCount[cacheIndex] == undefined) ? this.evalCount[cacheIndex]=1 : this.evalCount[cacheIndex]++;
		} catch (e) {
			// ignore exceptions
		}
		
		if ( this.optType < 0 ) {
			this.r[i] = -result;			// minimization: optimize for -result since RankVertices() assumes maximization
		} else if ( this.optType > 0 ) {
			this.r[i] = result;				// maximization: this is the default for RankVertices()
		} else {
			this.r[i] = -Math.abs(result);	// zero: optimize for -abs(result) since RankVertices() assumes maximization
		}
	}
	
	var ch;
	ch = "   ";
	switch (i) {
		case this.B: ch = "B  "; break;
		case this.N: ch = "N  "; break;
		case this.W: ch = "W  "; break;
		case this.R: ch = "R  "; break;
		case this.E: ch = "E  "; break;
		case this.Cr: ch = "Cr "; break;
		case this.Cw: ch = "Cw "; break;
	}
	
	// Display status
	if ( bShowStatus ) {
		var msg = "Working on " + this.name + "...\r";
		msg += "  Best Result: [";
		for ( d=0; d<this.nDim; d++ ) {
			msg += " " + this.bestSettings[d];
		}
		msg += " ]=" + Math.abs(-Math.round(this.bestResult*10000)/10000) + "\r";
		msg += "  Current Attempt: ["
			for ( d=0; d<this.nDim; d++ ) {
				msg += " " + this.uv[i][d];
			}
		msg += " ]=" + Math.abs(Math.round(this.r[i]*10000)/10000) + "\r";
		display.DisplayStatus(msg);
	}
}


// -------------------------------------------------------
// EvaluateResults()
// -------------------------------------------------------
Simplex.prototype['EvaluateResults'] = function (goal) {
	display.PostResults(this.name+":", (Math.abs(this.bestResult)<goal), Math.abs(this.bestResult));
	
	return Math.abs(this.bestResult);
}


// -------------------------------------------------------
// Initialize()
// Constructs an initial tilted simplex based on provided startingPoint
// which must be an array of nDim user-space coordinates.
// -------------------------------------------------------
Simplex.prototype['Initialize'] = function (startingPoint) {
	var d, i;
	var step, p, q;
	
	this.evalCount = new Array();	// vertex evaluation count for degenerate detection
	this.W = -1;					// index of Wastebasket vertex
	this.N = -1;					// index of Next-worst vertex
	this.B = -1;					// index of Best vertex
	
	for ( d=0; d<this.nDim; d++ ) {
		this.uv[0][d] = startingPoint[d];
	}
	this.NormalizeVertex(0);
	
	step = new Array();
	for ( d=0; d<this.nDim; d++ ) {
		step[d] = (this.v[0][d] > 0.75) ? (-0.25) : (0.25);
	}
	p = (Math.sqrt(this.nDim+1) + this.nDim - 1) / (this.nDim * Math.sqrt(2));
	q = (Math.sqrt(this.nDim+1) - 1) / (this.nDim * Math.sqrt(2));
	
	// Build the tilted simplex
	for ( i=1; i<this.nDim+1; i++ ) {
		for ( d=0; d<this.nDim; d++ ) {
			this.v[i][d] = this.v[0][d] + step[d] * ( (i-1 == d) ? (p) : (q) );
		}
	}
	
	// Evaluate the initial simplex
	for ( i=0; i<this.nDim+1; i++ ) {
		this.EvaluateVertex(i);
	}
	
	this.W = -1;		// index of Wastebasket vertex
	this.N = -1;		// index of Next-worst vertex
	this.B = -1;		// index of Best vertex
}

// -------------------------------------------------------
// Step()
// Perform an optimization step.  This may result in a move,
// expansion, or contraction of the simplex.
//
// Note:  If R is out of bounds, Cr is evaluated.  If Cr
// is in bounds, a Cr contraction is done rather than the
// default Cw contraction.  This is more in keeping with
// the spirit of the algorithm.
// -------------------------------------------------------
Simplex.prototype['Step'] = function () {
	var d;
	var useVertex;
	
	// step 1
	this.RankVertices();
	this.ComputeReflections();
	
	// step 2
	this.EvaluateVertex(this.R);
	if ( this.r[this.N] <= this.r[this.R] && this.r[this.R] <= this.r[this.B] ) {
		// We're moving in a mediocre direction.  Just accept it.
		// N <= R <= B    use simplex B..NR and go to step 3
		useVertex = this.R;
	} else if ( this.r[this.R] > this.r[this.B] ) {
		// We're moving in a desirable direction.  Try expanding.
		// R > B    calculate and evaluate E
		this.EvaluateVertex(this.E);
		if ( this.r[this.E] >= this.r[this.B] ) {
			// Expansion worked.
			// E >= B   use simplex B..NE and go to step 3
			useVertex = this.E;
		} else {
			// Expansion failed.
			// E < B   use simplex B..NR and go to step 3
			useVertex = this.R;
		}
	} else {
		// We're moving in an undesireable direction.  Try contracting.
		// R < N
		if ( this.r[this.R] == Number.NEGATIVE_INFINITY ) {
			// Boundary Violation
			this.EvaluateVertex(this.Cr); 	// evaluate Cr
			if ( this.r[this.Cr] != Number.NEGATIVE_INFINITY ) {
				// in bounds!   use B..NCr and go to step 3
				useVertex = this.Cr;
			} else {
				// left with a Cw retraction.   use B..NCw and go to step 3
				this.EvaluateVertex(this.Cw);
				useVertex = this.Cw;
			}
		} else {
			// No boundary violation
			if ( this.r[this.R] >= this.r[this.W] ) {
				// R >= W   calculate and evaluate Cr.  use simplex B..NCr and go to step 3
				this.EvaluateVertex(this.Cr);
				useVertex = this.Cr;
			} else {
				// R < W   calculate and evaluate Cw.  use simplex B..NCw and go to step 3
				this.EvaluateVertex(this.Cw);
				useVertex = this.Cw;
			}
		}
	}
	
	// step 3
	// discard W
	try {
		for ( d=0; d<this.nDim; d++ ) {	// for each dimension...
			this.v[this.W][d] = this.v[useVertex][d];
		}
	} catch (e) {
		alert(e + "\n useVertex=" + useVertex);
		alert("R=" + this.R);
		alert("E=" + this.E);
		alert("Cr=" + this.Cr);
		alert("Cw=" + this.Cw);
	}
	this.r[this.W] = this.r[useVertex];
	// transfer N to new W slot
	this.W = this.N;
	
	this.nSteps++;
	
	var ch, i, cacheIndex;
	for ( i=0; i<this.nDim + 1; i++ ) {
		ch = "   ";
		switch (i) {
			case this.B: ch = "B  "; break;
			case this.W: ch = "W  "; break;
		}
		this.DenormalizeVertex(i);
		cacheIndex = "";
		for ( d=0; d<this.nDim; d++ ) {
			cacheIndex += this.uv[i][d] + " ";
		}
	}
}

// -------------------------------------------------------
// GetElapsedMinutes()
// -------------------------------------------------------
Simplex.prototype['GetElapsedMinutes'] = function () {
	var now = new Date();
	
	return Math.round( ((now - this.startTime) / 1000 / 60) * 100) / 100;
}

// -------------------------------------------------------
// IsGoalReached()
// -------------------------------------------------------
Simplex.prototype['IsGoalReached'] = function () {
	var bGoalReached = 0;
	
	if ( this.optType < 0 ) {
		bGoalReached = (Math.abs(this.bestResult) < this.goal) ? 1 : 0;    // TODO: fix for negative goals
	} else if ( this.optType > 0 ) {
		bGoalReached = (this.bestResult > this.goal) ? 1 : 0;
	} else {
		bGoalReached = (Math.abs(this.bestResult) < Math.abs(this.goal)) ? 1 : 0;
	}
	
	return bGoalReached;
}

// -------------------------------------------------------
// IsDegenerate()
// Determines if simplex had degenerated to a single point
// or is oscillating
// -------------------------------------------------------
Simplex.prototype['IsDegenerate'] = function () {
	var i, d;
	var cacheIndex;
	var degenerate = 0;
	
	for ( i=0; i<this.nDim+1; i++ ) {
		this.DenormalizeVertex[i];
		cacheIndex = "";
		for ( d=0; d<this.nDim; d++ ) {
			cacheIndex += this.uv[i][d] + " ";
		}
		
		if ( this.evalCount[cacheIndex] > this.nDim*10 ) {
			degenerate = 1;
		}
	}
	
	return degenerate;
}

// -------------------------------------------------------
// Optimize()
// Begins optimization process.  Repeatedly calls Step()
// to do the optimization until one of the exit criteria
// defined by DefineOptimization() is reached.
// startingPoint must be an array of nDim user-space
// coordinates
// -------------------------------------------------------
Simplex.prototype['Optimize'] = function (startingPoint) {
	var d;
	var reason;
	
	display.DisplayStatus("");
	
	this.startTime = new Date();
	this.nSteps = 0;
	this.nRestarts = 0;
	this.bestResult = Number.NEGATIVE_INFINITY;
	
	// Evaluate starting point.
	for ( d=0; d<this.nDim; d++ ) {
		this.uv[0][d] = startingPoint[d];
		this.bestSettings[d] = startingPoint[d];
	}
	this.NormalizeVertex(0);
	this.EvaluateVertex(0);
	this.bestResult = this.r[0];
	
	// If goal reached on starting point, exit immediately.
	if ( this.IsGoalReached()  ) {
		reason = "Finished (Goal reached).";
	} else {
		
		try {
			this.Initialize(startingPoint);
		} catch (e) {
			alert(e);
		}
		
		// Optimize
		while ( 1 ) {
			if ( this.IsDegenerate() && this.B >= 0 ) {
				// simplex has degenerated to a point or is oscillating
				// restart using random points (best point is still retained)
				this.nRestarts++;
				
				if ( this.nRestarts <= this.restartLimit ) {
					for ( d=0; d<this.nDim; d++ ) {
						this.v[0][d] = Math.random();
					}
					this.DenormalizeVertex(0);
					this.Initialize(this.uv[0]);
				}
			}
			
			if ( this.IsGoalReached() && this.nRestarts > 0 ) {
				// goal reached and simplex converged at least once
				reason = "Finished (Goal reached).";
				break;
			}
			
			if ( this.GetElapsedMinutes() > this.timeLimit ) {
				// ran too long
				reason = "Terminated (Ran too long).";
				break;
			}
			
			if ( this.nSteps > this.stepLimit ) {
				// too many steps
				reason = "Terminated (Too many steps).";
				break;
			}
			
			if ( this.nRestarts > this.restartLimit ) {
				// too many restarts
				reason = "Terminated (Too many restarts).";
				break;
			}
			
			this.Step();
		}
	}
	
	for ( d=0; d<this.nDim; d++ ) {
		acr[this.fnSet[d]](this.bestSettings[d])
	}
	
	return this.bestSettings;
}
